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Abstract. Brenier and Grenier [SIAM J. Numer. Anal., 1998] proved that sticky particle dy¬ 
namics with a large number of particles allow to approximate the entropy solution to scalar 
one-dimensional conservation laws with monotonic initial data. In [arXiv: 1501.01498], we intro¬ 
duced a multitype version of this dynamics and proved that the associated empirical cumulative 
distribution functions converge to the viscosity solution, in the sense of Bianchini and Bres- 
san [Ann. of Math. (2), 2005], of one-dimensional diagonal hyperbolic systems with monotonic 
initial data of arbitrary finite variation. In the present paper, we analyse the error of this 
approximation procedure, by splitting it into the discretisation error of the initial data and the 
non-entropicity error induced by the evolution of the particle system. We prove that the error 
at time t is bounded from above by a term of order (1 + t)/n, where n denotes the number of 
particles, and give an example showing that this rate is optimal. We last analyse the additional 
error introduced when replacing the multitype sticky particle dynamics by an iterative scheme 
based on the typewise sticky particle dynamics, and illustrate the convergence of this scheme by 
numerical simulations. 


1. Introduction 

Systems of sticky particles have been known to reproduce the phenomenological behaviour of 
one-dimensional conservation laws in various physical contexts, in particular in astrophysics or in 
the study of gas dynamics [18, 16]. In such systems, finitely many particles evolve on the real line 
at constant velocity and stick together at collisions, with preservation of mass and momentum but 
dissipation of energy. The relation between these discrete systems and the equations of continuum 
physics was formalised by Brenier and Grenier [6], who showed that sticky particle dynamics with 
a large number of particles allow to approximate the entropy solution to scalar one-dimensional 
conservation laws with monotonic initial data. We also refer to Bouchut [4], Grenier [11], and E, 
Rykov and Sinai [8] for previous results in this direction. Based on this idea, we recently introduced 
a multitype sticky particle dynamics [13] in order to approximate the viscosity solution, in the sense 
of Bianchini and Bressan [2] , of one-dimensional diagonal hyperbolic systems with monotonic initial 
data of arbitrary finite variation. 

These sticky particle dynamics provide natural numerical schemes for the corresponding solu¬ 
tions to scalar conservation laws or diagonal hyperbolic systems. It is thus of interest to control 
the approximation error due to this procedure. This is the purpose of this article. We shall rely on 
the remark that sticky particle dynamics generically induce exact weak solutions to the considered 
equation, but for discrete initial data. Besides, these weak solutions need not satisfy Kruzkov’s 
entropy or Bianchini-Bressan’s viscosity condition. This leads us to split the total approximation 
error into a discretisation error of the initial data, and a non-entropicity error induced by the 
evolution of the particle system. 

The discretisation error of the initial data is addressed in Section 2. In particular, if the initial 
conditions have a compactly supported distributional derivative, this error in distance for n 
particles is proved to be bounded from above by a term of order 1/n. The error due to the 
evolution of the particle system is studied in Section 3 for the case of scalar conservation laws and 
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in Section 4 for the case of diagonal hyperbolic systems. In both cases, this error at time t > 0 
is proved to be bounded from above by a term of order tjn. This leads to a global convergence 
rate of order (1 + t)/n in the number n of particles. The precise statements for scalar conservation 
laws and diagonal hyperbolic systems are respectively given in Theorems 3.1 and 4.4, which are 
the main results of this paper. These results are finally illustrated with numerical simulations in 
Section 5. We emphasise the fact that the sticky particle approach is essentially restricted to the 
one-dimensional case and mention that (non-)existence and (non-)uniqueness issues related to its 
multidimensional generalisation were pointed out by Bressan and Nguyen [7]. 

The remainder of this introduction is dedicated to a detailed presentation of the sticky particle 
dynamics (SPD) and the multitype sticky particle dynamics (MSPD). 

1.1. SPD and scalar conservation laws. This subsection is dedicated to the introduction of 
the SPD, which allows to approximate the entropy solution to scalar conservation laws in one space 
dimension. 


1.1.1. Scalar conservation laws. Let us consider the scalar conservation law 


( 1 . 1 ) 


dfU -1- dx (A(r)) = 0, t > 0, X G M, 
u{0,x) = uo{x), 


for a nonconstant, monotonic and bounded initial condition uq. Up to an affine transform of 
the flux function A, one can assume that uq is the cumulative distribution function (CDF) of a 
probability measure m on the real line, which we denote uq = H ^ m where H{x) = l{a;>o} is the 
Heaviside function. The space of probability measures on the real line is denoted by P(M). Then 
A only needs to be defined on the interval [0,1], and it shall be assumed to have the following 
regularity. 

(C) The function A is of class on [0,1]. 

Under Assumption (C), we denote A = A' and Lq = sup^^[o,i] 

The following existence and uniqueness result follows from Kruzkov’s theorem, see [15, Theo¬ 
rem 2.3.5 and Proposition 2.3.6, pp. 36-37]. 


Theorem 1.1. Let A : [0,1] R satisfying Assumption (C) and uq = * m for m G P(M). 

There exists a unique weak solution u : [0,-hcxo) x M ^ [0,1] to the scalar conservation law (1.1) 
satisfying the entropy condition that, for all c G [0,1], 

dt\u -c\-\-dx (sgn(R - c)(A(u) - A(c))) < 0 

in the distributional sense, where sgr].{v) := l{a;>o} — l{cc<o}* 

In addition, it satisfies the following properties: 

(i) preservation of total variation: for all t > 0, u{t,') coincides dx-almost everywhere with 
the CDF of a probability measure rut on the real line; 

(a) finite speed of propagation: if uo{a) = 0, then u{t,a — tLc) = 0 for all t > 0, and if 
uo{b) = 1, then u{t, b -1- tLc) = 1 for all t > 0; 

(Hi) stability: if u and v refer to the entropy solutions to the scalar conservation law with 
respective initial data uq and vq, then for all t >0, 

') ~ '^(^5 OliLqR) < ll'^o - '^oIIlI(R)- 

1.1.2. Sticky Particle Dynamics. For n > 1, we denote by D^ the polyhedron of RF defined by 

:= {x = (xi,..., Xn) : < • • • < 

Let X G Dn be a vector of initial positions and A = (Ai,..., A^) G a vector of initial velocities. 
Under the SPD, the particle with index k G {1,... ,n} has initial position x/e, initial velocity 
and mass 1/n. It evolves at constant velocity on the real line, up to the first collision with another 
particle. At collisions, the particles stick together and form a cluster: its mass is given by the 
number of colliding particles over n, and its velocity by the average of the pre-collisional velocities 
of the particles. More generally, when several clusters collide, they form a single cluster with 
conservation of total mass and momentum. 
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For all t > 0, the position of the k-th particle at time t > 0 is denoted by 0/c[A](x;t), and it 
is easy to check that the process (0[A](x; t))t>o defined by 0[A](x;t) := ((/)i[A](x; t),... ,0n[A](x;t)) 
induces a continuous flow in Dn- Its stability with respect to the initial configuration and the 
vector of initial velocities is detailed in Proposition 1.2 below. Before stating this result, let us 
define the normalised distance on by 

l|x-y||i := -yfcl- 

/c=l 

Proposition 1.2. [13, Proposition 3.1.9, (i)] Let x, y G and A,^ G For all 0 < s <t, 

, _ n 

(1.2) ||<(>[A](x;i) -(/)[^](y;t)||i < ||(/)[A](x; s) - (?^[/l](y; s)||i + - 

k=l 

1.1.3. Approximation of the sealar eonservation law. Let A satisfy Assumption (C). In order to 
approximate the entropy solution to the scalar conservation law (1.1), we specify a choice of initial 
velocities for the SPD by defining A G as 

pkin 

V/c G {1,..., n}, Xk ’= n / A(re)dre. 

J w={k—l)/n 

Given an initial configuration x G we define the empirical distribution of the SPD at time 
t > 0 by 

1 "" 
k=l 

and the associated empirical CDF by 

1 "" 

Un[A{t,x) ■.= H * = -X!^W.[A](x;t)<^}- 

k=l 

Given n > 1 and x G it is easily checked that the empirical CDF i4n[x] satisfies the 
properties (i), (ii) and (hi) of Theorem 1.1. The preservation of the total variation is obvious, 
and the finite speed of propagation is just the transcription of the fact that the modulus of the 
initial velocities Xk is bounded by Lc, uniformly with respect to n. Finally, the stability follows 
from (1.2) with jh = A. 

It can also be shown that is a weak solution to the scalar conservation law (1.1) with 

discrete initial condition [13, Proposition 4.2.1]. However for fixed n, it does not necessarily satisfy 
the entropy condition of Theorem 1.1. This property is recovered when taking the limit of an 
infinite number of particles, as is expressed by the following result, the proof of which is originally 
due to Brenier and Grenier [6], see also [12] and [13, Lemma 8.2.3] for appropriate generalisations. 

Theorem 1.3. Let A satisfy Assumption (C) and let m G P(M). Let (x(n))n>i be a sequenee of 
initial eonfigurations sueh that, for all n>l, x(n) G and the empirieal distribution 

1 "" 

Mo[x(n)] = 

^ k=l 

eonverges weakly to m. 

For all t > 0, the empirical distribution jat[x{n)] eonverges weakly to the probability measure 
rut G P(M) sueh that u(t,x) := H ^mt{x) is the unique entropy solution of the sealar conservation 
law (1.1) with initial condition uq = F[ ^ m. Equivalently, for all t > 0, the empirical CDF 
Un[x{n)]{t, •) converges dx-almost everywhere to u{t, •). 

Using Theorem 1.3 to pass to the limit n +oc in (1.2), the stability inequality of (hi) can be 
extended in order to take the dependence of the entropy solution on the flux function into account. 
This is done in the next proposition, which is of independent interest, and the proof of which is 
postponed to Appendix A. 
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Proposition 1.4. Let A, M : [0,1] ^ M satisfying Assumption (C), and be CDFs on the 

real line. Denote A := A' and /i := M', and call u and v the entropy solutions of the scalar 
conservation law with respective flux function A and M, and respective initial data uq and vq. 
Then, for all 0 < s <t, 

\\u{t, •) - v{t, OIIlMR) < ll«(s> •) - OIIlHi;) + (i - s) / \Hw) - iJ,{w)\dw. 

J w=0 

1.1.4. Rate of convergence. Given a sequence of initial configurations (x(n))^>i satisfying the 
assumptions of Theorem 1.3, our first purpose in this article is to estimate the error when ap¬ 
proximating u with R^[x(n)]. On account of the stability property stated in Theorem 1.1, a fairly 
natural distance to measure this error is the distance on M. Indeed, this stability property 
allows us to write, for all t > 0, 

\\u{t,-) -M„[x(n)](i, OIIlhm) 

(1.3) < ||■u(^, •) - ■Uoo[x(n)](i, OIIl^k) + ll'Woo[x(n)](t, ■) - ■u„[x(n)](i, OIIl^k) 

< \\uo - Moo,o[x(n)]||Li(E) + ||woo[x(n)](t, •) - w„[x(n)](i, •)||li(R), 

where we have introduced the entropy solution UoQ[x{n)] to the scalar conservation law (1.1) with 
discretised initial condition Roo,o[x(^)] := */io[x(n)]. The two terms in the right-hand side above 

are of a very different nature, and can be estimated separately: the first term corresponds to the 
discretisation error of the measure m, while the second term only measures the non-entropicity 
induced by the evolution of the particle system for a given initial condition Roo,o[x(n)]. 

The discretisation error of the measure m is addressed in Section 2. There, we use the fact that, 
given two probability measures m and m' on the real line, the distance between H ^ m and 
* m' is the Wasserstein distance of order 1 between m and m ', defined by 

(1.4) Wi(m,m') := inf / \x — x'\m{dxdx'), 

where the infimum runs over all the probability measures m G P(M^) such that, for all Borel sets 
A, A' c M, 

m(A X M) = m{A), m{R x A') = m'{A'). 

This is due to the fact that, on the real line, the measure 

m = U o ((iJ * to)-\ {H * , 

where U refers to the Lebesgue measure on [0,1], realises the infimum in (1.4). In this definition, 
the pseudo-inverse F~^ of a CDF is defined by 

(1.5) V'z; G (0,1), F~^{v) := inf{x G M : F{x) > v}. 


We deduce that 

(1.6) Wi(m,m')= f \{H ^ m)~^{v) — {H ^ m')~^{v)\dv = f \H ^ m{x) — H ^ m'{x)\dx, 

Jv=0 JxER 

whence \\H *m — H * m'||Li(R) = Wi(m, m'). 

As a consequence, the first term in the right-hand side of (1.3) rewrites 


\\uo - Woo,o[x(n)]||Li(R) = Wi(m,/u”), 


h" 




{jd " 


k=l 


Precise bounds on Wi(m,/i’^) in terms of n for the optimal discretisation of m are derived in 
Lemma 2.1. They depend heavily on the tail of m. In particular, an important remark to be done 
at this point is the following. Assume that m has an infinite first order moment. Then, since 
by (1.4) and the triangle inequality. 


Wi(m,m') > 


/ \x\m{dx) — / 
JxER Jx'ER 


\x'\m' {dx') 
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any approximation of m by a measure with finite first order moment necessarily satisfies 
Wi= + 00 . As a consequence, there is no purpose in trying to compute a rate of con¬ 
vergence in this case. Therefore, although our results hold true without any assumption on m, 
they only have a nontrivial content when m has a finite first order moment. 

The non-entropicity error is then addressed in Section 3, where given arbitrary n > 1 and 
X G an estimation is first derived on the distance between rt^[x](t, •) and i 4 oo[x](t, •) in 
Proposition 3.2. This result holds under the following strengthening of Assumption (C). 

(LC) The function A = A' is I/Lc-Lipschitz continuous. 

Combining the results of Section 2 with Proposition 3.2 yields complete rates of convergence of the 
SPD, as is stated in Theorem 3.1. 

1.2. MSPD and diagonal hyperbolic systems. This subsection is dedicated to the introduc¬ 
tion of the MSPD, which allows to approximate the semigroup solution to diagonal hyperbolic 
systems in one space dimension. We refer to [13] for more details. 

Let us fix an integer d > 2 , and consider the diagonal hyper- 

J dtu^ + \^{\x)dxU^ = 0 , t > 0 , X G M, 

=ul{x), 

for nonconstant, monotonic and bounded initial data itj,..., Once again, we shall assume that 
there exists m = (m^,..., m^) G P(M)^ such that, for all 7 G {1,..., d}, = H and look for 

solutions u = ..., u^) of (1.7) such that, for all t > 0, for all 7 G {1,..., d}, u^{t^ •) remains 

the CDF of a probability measure on the real line. The characteristic fields A^,...,A^ are 
therefore defined on [ 0 , 1 ]^, and we shall extend Assumptions (C) and (LC) as follows. 

(C) For all 7 G {1,..., d}, the function is continuous on [0,1]^. 

Under Assumption (C), we denote Lq = maxi<^<(^sup^^^o,!]^ 

(LC) There exists Llc ^ [0, + 00 ) such that 

d 

V 7 e {l,...,d}, Vu,v e [ 0 , 1 ]“^, |A^(u) - A^(v)| < Llc X] 

7 ' = 1 

We also require the system to be uniformly strictly hyperbolic, in the sense of the following as¬ 
sumption. 

(USH) There exists Tush ^ (0, + 00 ) such that 

V7 e 1}, inf A'>'(u) - sup A'>'+^(u) > Lush- 

uG[0,1]^ 

On account of the fact that the system (1.7) is written in a nonconservative form, both the notions 
of weak and entropy solution are not canonically defined. An appropriate notion of weak solution is 
introduced in [13, Definition 2.4.1], while a criterion for uniqueness is stated in [13, Definition 8.2.5] 
by adapting the notion of viscosity solution by Bianchini and Bressan [2], to which we shall refer 
as the semigroup solution to (1.7). These notions are used in Theorem 1.6 below. 

1.2.2. Typewise and Multitype Sticky Particle Dynamics. In order to approximate each coordinate 
of the solution to the system (1.7), we shall now introduce d systems of n particles evolving on 

the real line, each system being associated with a type 7 G {1,..., d}, so that the empirical CDF 
of the system of type 7 is supposed to approximate . The k-th particle of type 7 is referred to 
as the particle 7 : /c, and the set of such indices is denoted P^. A configuration is described by an 
element x = of the Cartesian product D^. The normalised L^ distance is defined on 

Dt by 


1 . 2 . 1 . Diagonal hyperbolic systems. 
bolic system 

(1.7) V7 e{i,...V}, 
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The first step towards the definition of the MSPD is the introduction of the Typewise Sticky 
Particle Dynamics (TSPD). Given an array of initial velocities A = and an initial 

configuration x G we denote by 


l>[A](x;d = (l>^[A](x;t))^,fegpd 

the process in obtained by letting the system of type 7 evolve according to the SPD with initial 
configuration G Dn and initial velocity vector = (A^,..., A^) G without 

any interaction with other systems. 

The MSPD is built from the TSPD as follows: 

(i) the particle j : k is initialised with an array of velocities corresponding to a discretisation 
of the function A^ at a point u G [0,1]^ recording the rank of the particle 7 : /c in each of 
the d systems, see ( 1 . 8 ) below; 

(ii) when clusters of particles of different types collide, they cross each other and the TSPD 
is restarted with initial velocities depending on the post-collisional rank of the particles in 
each system. 

Assumption (USH) essentially implies that whatever the arrangement of the particles, particles 
of lower type always have a larger velocity. This prescribes the post-collisional order to update 
the velocities of the TSPD, and ensures that clusters of different types drift away from each other 
immediately after a collision. For an initial configuration x G the array of initial velocities 
M^) = (7W)7 :keP^ is defined under Assumption (C) by 


( 1 . 8 ) 


— 


pk/n 

' / 

J w=(k- 


w=(k — l) / n 



,7-1 

^7:/c 


(x) 


re, uj 


7+1 

j-.k 



die. 


where u;(J^.^(x) denotes the (scaled) rank of the particle 7 : k within the system of type 7 ', formally 
defined by 


^j-.k 


(x) := 


n 2-^k' = l 
I 

n = l 




if 7 ' < 7 , 
if 7 ' > 7 , 


where particles of different types sharing the same location are counted according to the post- 
collisional rank imposed by Assumption (USH). 

The resulting dynamics in is the MSPD, denoted by ($(x;t))t>o- More details on its 
construction are given in [13, Subsection 3.2], where it is proved that it defines a continuous flow. 
Its stability with respect to the initial configuration is described by the next result, which forms 
the core of the article [13]. 


Proposition 1.5. [13, Theorem 2.5.2] Under Assumptions (LC) and (USH), there exists Ci G 
[1, -hoc) depending only on d and the ratio Tlc/Aush sueh that, for all n > 1, for all x, y G 
for all 0 < s <t, 

||$(x;f) - $(y;f)||i < £i||$(x;s) - $(y;s)||i. 


In contrast with Proposition 1.2, we note that no stability result with respect to the characteristic 
fields A^,..., A^ is available. 


1.2.3. Approximation of the diagonal hyperholie system. Similarly to the scalar case, we define the 
empirical distribution of the system of type 7 in the MSPD at time t > 0 as the probability measure 
on the real line given by 

n 

k=l 

and the vector of empirical CDFs Un[x] = (r;^[x], ..., r^[x]) by 

1 ^ 7 

k=l 
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It is proved in [13, Proposition 4.2.1] that, for all x G Un[x] is a weak solution to the 
system (1.7). Given a sequence of initial configurations (x(n))n>i such that, for all n > 1, x G 
and, for all 7 G {1,..., d}, the empirical distribution 

Aio[x(n)] = 

k=l 

converges weakly to some G P(M), one could by analogy with Theorem 1.3 expect Un[x(n)] to 
converge to a weak solution of the system (1.7) with initial data Uq=H^ rri^ ^ satisfying some 
specific entropy-like condition making it unique and physically meaningful. Although it is true that, 
up to extracting a subsequence, u^[x(n)] actually converges to a weak solution of the system [13, 
Theorem 2.4.5], following Bianchini and Bressan’s construction we were only able in [13] to identify 
the limit if it satisfies some semigroup and stability estimate with respect to the initial data. For 
this purpose, we introduced [13, Definition 2.6.4] a discretisation operator Xn • P(l^)^ ^ 
defined by = x(n) with 

.( 2 /c+l)/( 2 (n+l)) 

(1.9) Vy :/c G P^, x^(n) = (n + l) / {H ^ m^)~^{w)dw. 

J w=(2k—l) / ( 2 (n+l)) 

For all 7 G {l,...,d}, the empirical measure converges weakly to m^, see [13, 

Lemma 8.1.5], and we obtained the following convergence theorem for the MSPD, which follows 
from the more general statement of [13, Theorem 2.6.5]. 

Theorem 1 . 6 . Under Assumptions (LC) and (USH), let m = (m^,..., m^) G P(M)^. 

For all t >0, for a// 7 G {1,..., d}, the empirical distribution converges weakly to the 

probability measure G P(M) such that u = (i 4 ^,... defined by u^{t^x) = P * m]{x) is the 
unique semigroup solution to the diagonal hyperbolic system (1.7) with initial data Uq=H^ . 
Equivalently, for all t > 0, for all 7 G {!,..., d}, the empirical CDF rtj^[x^m](t, •) converges 
dx-almost everywhere to u^(t, •). 

Each coordinate of the semigroup solution u satisfies the preservation of total variation and 
finite speed of propagation properties stated in Theorem 1.1. The stability property writes as follows: 
for all m, m' G P(M)^; the semigroup solutions u and v to the system (1.7) with respective initial 
data uq; vq defined by Uq = H ^ m?, i;q = P * m!^ satisfy, for all 0 < s <t, 

d 

||u(i, •) - w{t, •)I|li(E)‘! := T. •)IIli(E) 

7 = 1 

d 

•) - IIIlHE) =; A||u(s, •) - v(s, •)IIl 1 (E)‘*- 

7=1 

The notion of semigroup solution is adapted from [2] and detailed in [13, Definition 8.2.5]. 

In Corollary 4.6, we shall generalise this result and prove the convergence of the empirical cumu¬ 
lative distribution functions of the MSPD to the semigroup solution for a large class of sequences 
of initial configurations (x(n))n>i approximating a vector of initial measures m = (m^,... ,m^) 
such that, for all 7 G {1,..., d}, mC has a finite first-order moment. 

1.2.4. Rate of convergence and approximation by the iterated TSPD. The second purpose of this 
paper is to supplement Theorem 1.6 with a rate of convergence of u^[x^m](t, •), or more generally 
Un[x(n)](t, •) for a sequence of initial configurations (x(n))n>i approximating m, to u(t, •). As in 
the scalar case, we split the approximation error by writing 

||u(t, ■) - u„[x(n)](i, •)||Li(M)<i 

< \\u{t, ■) - Uoo[x(n)](t, OllLnE)-* + l|uoo[x(n)](t, •) - u„[x(n)](i, •)IIli(E)<‘ 

< >Ci||uo - Uoo,o[x(n)]||Li(R)c! + |uoo[x(n)](i, •) - u„[x(n)](t, •)||li(e)‘*) 

where Uoo[x(n)] is the semigroup solution, in the sense of Theorem 1.6, of the system (1.7) with 
initial condition Uoo,o[x(n)] defined by, for all 7 e *Aio[x(n)]. In the 

last line above, we used the stability estimate of Theorem 1.6 A on semigroup solutions. 
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The estimation of the discretisation error of the initial data follows from the analysis of the 
scalar case carried out in Section 2 . We shall use the distance on P(]R)^ defined by, for all 

m = (m^,..., m^), m' = (m'^,..., m'^) G P(M)^, 

d d 

(1.10) W^''^(m,m') = y]Wi(m^,m'^) = y] ||iJ * * m''>'||Li(E) = ||u - v||Li(R)<i, 

7=1 7=1 

with u := (i^ * • • • , * m^), v := (i^ * m'^, • • • , * m'^). 

The analysis of the error due to the evolution of the MSPD turns out to be more delicate than in 
the scalar case, and we shall resort to an approximation of this dynamics by the following scheme. 
Fix a time step A > 0 and define the process (4 >a(x; t))t>o in letting, on each time interval 

[{L — 1 )A,LA], 1/ > 1 , the particles evolve according to the TSPD with initial velocity vector 
A(4 >a(x; {L — 1 )A)). This amounts to neglecting the collisions between clusters of different types 
on [{L — 1 ) A, LA] and only updating the velocities with respect to the new ordering of the system 
at the end of each such interval. From a computational point of view, this approximated dynamics 
is expected to be easier to simulate than the MSPD, as one does not have to keep track of all the 
collisions. The approximation error induced by this scheme is addressed in Section 4, where we 
also use it as a theoretical tool to obtain rates of convergence in Theorem 4.4. We finally present 
a numerical implementation of this scheme in Section 5. 


2 . Wi discretisation error of initial conditions 

In this section, we want to approximate the probability measure m on the real line, with CDF 
F = * m, by the empirical measure 




n 

n 


k=l 


of a vector x = (xi,... ,Xn) of n deterministic points on the real line, with xi < • • • < On 
account of ( 1 . 6 ), one has 

n nF-^(k/n) 

Wi(m,/i") = V/ \y-Xk\m{dy), 

k=l Jy=F-^iik-l)/n) 

where we recall the definition (1.5) of the pseudo-inverse F~^ and complement it with the conven- 
tion F-i(O) = in7g(o,i) and = sup„g(o,i) F~'^{v). 

The choice 

( 2 . 1 ) 


Xk = F 


-1 


2k-1 
2n 


the median of the image of the uniform law on [{k — 1)/r, k/n] by F for all k G {1,... ,r}, 
minimises this Wasserstein distance. Let us now analyse this optimal choice. First, 

^ nk/n 

Wi(m,//") = V / (F-\u) - F-\u-l/{2n)))du 

J U=(^ .. 


k=l '^u=(2k — l)/(2n) 


< 


When 


f (F- 1 (m)-F- 1 (m- l/( 2 n)))dw. 

iu=l/(2n) 

F~^{u)\du < +OC, 


/ 

J u=0 


that is to say m has a finite first order moment, one deduces that 

nl /*l/(2n) 

Wi(m,/i’^) < / F~^{u)du— / F~^{u)<lu^ 

Ju=l — l/{2n) J u=0 

where the right-hand side goes to 0 when n +oo by Lebesgue’s theorem. With the discussion 
of the case where the first order moment of m is infinite at the end of §1.1.4, we deduce the first 
statement in the next lemma. 
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Lemma 2.1. Let nL he defined as in the discussion above. The behaviour 0 /Wi(m,/i’^) depends 
on the tail of m as is described below. 

(i) One has lim^^+oo Wi(m ,= 0 if and only if m has a finite first order moment. 

(a) For all n > 1, 

^ [ ^yF{x){l - F{x))dx, 

JxeR 

where the right-hand side may be infinite. 

(Hi) If there exist —00 < a < b < +00 such that m([a, b]) = 1, then for all n > 1, 

In 

(iv) If m{dx) = ^{xe[a,h\}f{x)^x / positive on [a, 6] where —00 < a <b < + 00 ^ then 

lim nWi ^ 

n ^+00 4 


Proof. The assertion (ii) follows from the comparison of with the expected Wasser- 

stein distance between m and the empirical measure ^ Ylk=i of independent random variables 
Xi,..., Xn with identical distribution m, a detailed study of which was carried out by Bobkov and 
Ledoux [3]. By the optimality of the choice of we first have 


Wi< E 


Wi 



1 

n 



According to the proof of [3, Theorem 3.2], the right-hand side is not greater than 


£1/2 


= 

'(/ 


\ n ^' / 

L \ k=i / J 


\J 


^^{Xk<x} 


k=l 


dx 


f 1 ^ 


< f 

^ [ ^Fix)il-Fix))dx. 

/n JxeR 


{Xk<x} 


dx 


If m([a, b]) = 1, then the mass of the measure u on (0,1) such that z^((0, u)) = F ^{u) — F ^(0+) 
for all 14 G (0,1) is smaller than b — a. Moreover, 


Wi 


it p 

k = l’^^ 
n p 

=s/ 

u -\ ’J V' 


kin 

Ju={2k-l)/{2n) Jv=' 


/ '^{u-l/( 2 n)<v<u}v{dv)du 

Jv=0 


< 


Jve[{k-l)/n,k/n) 

K( 0 d)) _ 

2n 


mm 1 v — - - — — v] h>(dv) 

n n 


Let us finally assume that m{dx) = ^{xe[a,h\}f{x)^x with / positive on [a, 6]. Since m is the 
image of the Lebesgue measure on [0,1] by one has 


F-\u) 


-\u)-F-\u-ll{2n))= [ 

Jx=F-^ 


m{dx) 


f 

J v=u 


di’ 


(u—l/(2n)) /(^) Jv=u—l/(2n) f (^)) 

for all u G (l/(2n), 1). The continuity of translations in ensures that 

1 


lim [ 

Ju=0 


nt{ue(i/( 2 n),i)}{F ^{u)-F - l/(2 n)))- 


2f{F-\u)) 


du = 0. 
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2 ^ J«=(2fc-l)/(2n) f{P W) 

- /_o ”^{«e(i/(2n),i)} - l/(2n))) - 


it is enough to check that 


lim > 


du 1 


2 y„=o /(J’-'N)' 

This follows from the weak convergence of 2 l{ue((2/c-i)/(2n),/c/n)}dR to l{ne(o,i)}dR and the 

density of continuous and bounded functions in L^(]R). □ 

Remark 2.2. Since 


the condition 


is equivalent to 


Similarly, the condition 


' \F ^(R)|dR= / \x\m{dx) = / {{1 — F{x))F{—x))dx, 

u=0 JxEM. J x=0 


\F ^{u)\du < +00 


F{x){l — F{x))dx < + 00 . 


^F{x){l — F{x))dx < +00 


is an assumption on the decay of the tails of m. Since F{x){l — F{x)) < ^/F{x){l — F{x)), it 
implies that 


On the other hand, if for some a > 2, 


\x\m{dx) < + 00 . 


\x\'^m{dx) < + 00 , 


(1 + ^)F{x){l — F{x))dx < Too 


and, by the Cauchy-Schwarz inequality. 


( / ^/F{x){l — F{x))dx \ < / (1 + ) ^dx / (1 + )F(x)(l — F(x))dx < Too. 

\JxeM / JxeM. JxeM. 

We finally discuss some other cases where Wi(m,/i^) can be estimated. If m has a positive 
density / on ( 0 ,oo), one has 


Wi (m,/i-) = 


k=l v={k-l)/n 
^ rl-l/{2n) 


. f k — 1 k \ di; 

mm V -,- V ] ^ 

V n n ) f{F-^{v)) 


1 .1 (l-v)dv 

- ^ X=o f{F-\v)) + ii_a/(2„) f{F-\v)) 


F-^ 1 


F-\0-^ 


(1 — v)dt> 

= l-l/(2n) f{F-^v)y 
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On the other hand, if u 1/f{F is nondecreasing on {u, 1 ) with u G (0,1), then 


Wi 


nk/n 

= J 2 

7.-1 Jv=(k- 


k — 1 k 
mm I V -,- V 


k=l ^v=ik-l)/n 
^ y^(4/c —l)/(4n) 


>^y r 

4n^A= 


n n 
dv 


dv 


f{F-Hv)) 


^n=(4/c—3)/(4n) f (^)) 


s / 

_, r /A 'J V 


(4/c-l)/(4n) 


dt> 


^ 1 


>-[ 

8 n y„ 


l-l/4n 


di; 


n+l/4n f{F-Hv)) 


1 


8n 




An 




An 


We apply these estimates on the following two examples. 
• For F{x) = “ x~^) with ct > 1, one has 

= (1 - 

Since 


f 

Jv=l- 


(1 — v)dv 


= 0{n 


-l+l/a 


.=l-l/(2n) fiF-^v)) 

and for c > 0 , F~^{1 — c/n) = 0{rdl^), we conclude that Wi(m,/r’^) = 
• For F{x) = ~ exp(—with ct > 0, one has 

F-i(w) = (-log(l-«))!/“, 

and we conclude that Wi= 0 ((logn)^/^/n). 


3. Rate of convergence of the SPD to scalar conservation laws 

Under the assumptions of Theorem 1.3, the purpose of this section is to estimate the distance 
between the entropy solution u of the scalar conservation law ( 1 . 1 ) with initial condition lig = 
for some m G P(M), and the empirical CDF ix^[x(n)] associated with the SPD started at some 
configuration x(n) = (xi(n),... ,Xn(n)) G Dn^ over finite time horizons. 

Theorem 3.1. Let A satisfy Assumption (LC). Then for all n > 1, for all x(n) G Dn, 

\/t > 0, \\u{t, •) - u„[x(n)](i, •)||li(e) < Wi(m,/Lio[x(n)]) + 

When m([a, b]) = 1 with — oo < a < b < +oo and x(n) = (xi(n ),..., Xn(n)) is given by 

f2k — 1\ 

(3.1) Xkin) = j , ke{l,...,n}, 

then for all n > 1, 

Vi > 0, \\u{t, ■) - Un[x{n)]{t, OIIlhm) < -— 

Proof Following the approach described in the introduction of the article, we first write, for all 
t > 0, 

,^ 2 ) -Mn[x(n)](i,.)||Li(M) < ||«(i,-)-Moo[x(n)](i,-)||Li(M) 

+ ||woo[x(n)](i, •) - Un[x{n)]{t, •)||li(E), 

where iXoo[x(n)] refers to the entropy solution of the scalar conservation law (1.1) with discretised 
initial condition rtoo,o[x(n)] = H ^ /io[x(n)]. The stability property of Theorem 1.1 for the scalar 
conservation law (1.1) yields 

||M(i, ■)-Moo[x(n)](i, ■)||li(b) < Iko - Moo,o[x(n)]||Li(R) = Wi(to,/ 7o[x(n)]), 
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which, according to (hi) in Lemma 2.1, is smaller than {b — a)/(2n) if m([a, 6 ]) = 1 and x(n) is 
given by (3.1). The following proposition allows to control the second term in the right-hand side 
of (3.2). □ 


Proposition 3.2. Under Assumption (LC), let n > 1 and let x G D^. For all t >0, 


l|Moo[x](t, 




tLi,c 

n 


Proof. Instead of comparing [x] directly with the entropy solution [x] , we first compare it with 
the empirical CDF of the SPD with 2n particles, started in the configuration where the positions 
in X are duplicated, and then iterate this comparison. More precisely, for all n > 1, let us define 
the operator 

Dn D2n 

X X 

by X 2 k-i = X 2 k = Xk for all /c G {1,... ,n}, and denote by x"^ G D 2 mn the m-th iteration of this 
operator. We first prove that, for all n > 1, for all x G 

(3.3) \/t > 0, ||Mn[x](t, ■) - M2n[x](t, •)||li(R) < ' 

In this purpose, we interpret the function Un[^] as the empirical CDF of the SPD with 2n particles, 
with initial configuration x, and with modified velocity coefficients Ai,..., X 2 n defined by 


A 2 /C -1 = ^2k =n 

J u 

for all /c G {1,..., n}. Then (1.2) yields 


•kin 
w={k — l)/n 


X{w)dw, 


2n 


lkn[x](t, •) - W 2 n[x](i, OIIli - fP, 

k=l 

n 


Xf. - 2n 


•kjn 


rk/{2n) 

/ A(R;)dR; 

J w={k — l)/{2n) 


w={k — l)/n 
'kjn 


Mk-^)tn 

X{w)dw — 2 / A(R;)dR; 

J w={k — l)/n 

nkjn nkjn ^ 

/ X{w)dw — 2 / A(R;)dR; > 

Jw={k—l)/n Jw={k—^)/n J 


and Assumption (LC) allows to bound the right-hand side above by tI/Lc/(2n), whence (3.3). 
We now fix n > 1 and use (3.3) to write that, for all M > 1 , 


M 


||M„[x](i, ■) -Moo[x](t, OIIlHR) ^ E ll'“2—inlx™ (i,[x™] (i, •) ||l1( 

m=l 

+ \\U2Mn[x^]{t, •) - Uoo[x](t, OIIl^R) 


M 

- E ^i^ + ll“ 2 M„[x^](i,-)-Moo[x](i,-)l|Li(R)- 


m=l 


To complete the proof, we have to check that 
(3.4) -- 


lim \\u 2 M„[x^]{t, ■) - Moo[x](i, •)IIli(r) = 0 . 


To this aim, we note that, for all M > 1, the empirical measure associated with x^ does 

not depend on M and coincides with /ioH- As a consequence, by Theorem 1.3, R 2 ^n[^^](C') 
converges to Roo[x](t, •), dx-almost everywhere. Besides, the initial measure /ioH has compact 
support {xi,...,x^}. Let —oc < a < xi and Xn < b < +oo, so that Roo,oH(<^) = 0 aiid 
'^ 00,0 H(^) = 1- Using the finite speed of propagation for both R 2 m^[x^] and RooH, we deduce 
that 

\fx ^ [a-tLc,b^tLc], U 2 M^[x^]{t,x) -Roo[x](t,x) = 0. 

Therefore, R 2 M^[x^](t, •) — Roo[x](t, •) has a compact support which does not depend on M. As a 
consequence, (3.4) follows from Lebesgue’s theorem and the proof is completed. □ 
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Remark 3.3. If the flux function A is concave, then itn[x] is actually the entropy solution to the 
scalar conservation law (1.1) with discrete initial condition [5, Lemma 3.3], so that 


||■Woo[x](t •) - Un[x]{t, •)||l1(R) = 0, 


for all n > 1, X G and t > 0. Note that this remark holds true even without Assumption (LC). 

On the contrary, one can construct examples of convex flux functions for which the bound of 
Proposition 3.2 provides the right order of magnitude for ||iXoo[x](t, •) — 'Un[x](t, •)||li(r). Let us 
indeed consider the Burgers equation A{u) = ix^/2 with the initial condition m = the Dirac 
mass at 0. Fix n > 1 and let x = (0,..., 0) G D^. On the one hand, the entropy solution iXoo[x] 
writes 


Moo [x] (i, X) 


0 if X < 0, 

< I if 0 < X < t, 
1 if X > t. 


On the other hand, the particles drift away from each other in the SPD, so that 


Vfc e (pk[>^]{x;t) =t 


fc - 1/2 


therefore 

Un[x]{t,x) = < 

As a consequence. 


0 ifx<^, 

- if < X < fh+m with k E {1 ,..., n — 1}, 

— n 


||Woo[x](t, •) - Wn[x](t, OIIlHR) = 

which is of the same order as the bound given in Proposition 3.2. 


4. Rate of convergence of the MSPD to diagonal hyperbolic systems 

4.1. Approximating the MPSD through the iterated TSPD. Given A > 0, we define the 
operators 4 >a and 4 >a on by 

$a(x) := i'(x; A), ^a(x) := 4-[A(x)](x; A), 

where we recall that 4>(x; t) refers to the MSPD while 4>[A(x)] (x; t) is the TSPD with initial velocity 
vector A(x) given by (1.8). 

By the flow property of the MSPD, the L-th iteration 4>£(x) is nothing but 4>(x;I/A). On the 
other hand, the L-th iteration 4 >a(x) is easier to compute as one does not have to take interactions 
between particles of different types into account. The purpose of this subsection is to estimate the 
error of the approximation of the MSPD by this iterated TSPD scheme. To this aim, given x G 
we first extend the iterated TSPD into a continuous process (4>a(x; t))t>o in by interpolating 
between the points of the uniform grid with step A thanks to the TSPD. More precisely, for all 
L > 1, for alH G [{L — 1)A, LA], we define 

$A(x;f) :=$*_(z,-i)A($i-'(x)) = $[A($i-'(x))]($i-i(x);f-(L-l)A). 

Our purpose is to prove the following result. 

Proposition 4.1. Under Assumptions (LC) and (USH), there exists C > 0 that depends neither 
on n>l nor on A such that, for all x G D^, 

(4.1) sup ||4>(x;t) - 4>A(x;t)||i < CA. 

t>o 

The value of C is given in (4.6) below. Of course, for t restricted to the uniform grid (I/A)i,>o, 
the estimate (4.1) implies 

sup||i>i(x)-$i(x)||i ^CA. 

L>0 
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We begin by proving this restricted estimation before deducing (4.1). We first use the flow property 
of both the TSPD and the MSPD [13, Proposition 3.2.8] to write 

II^a(x) - $i(x)||i < ||$i-'($A - $A)$i-'(x)||i 

(4.2) 

^=1 

where the last inequality follows from the discrete stability estimate of Proposition 1.5. 

The next lemma allows to estimate each term of the sum in the right-hand side. We recall 
from [13] that, for all y G the set R(y) contains the pairs of particles {a : ^ : j) ^ (^n)^ such 

that a < f3 and yf < Vj ^ so that in the MSPD started at y, these particles undergo a collision at 
a finite and positive time We also define the first collision time in the MSPD started 

at y by 

£(y) := min{T™“^y (y), {a : i, p : j) € R(y)} € ( 0 , +oo], 
where we take the convention that t*(y) = +oc if R(y) is empty. 

Lemma 4.2. Under the assumptions of Proposition 4-1, for all y G and all A > 0, 

II A(y) ~A(y)||i _ ^2 E 

(a:i,/3:i)eR(y) 


Proof of Lemma 4-2. Let to := 0 < ti < • • • < ti? < A =: tRj^i refer to the successive instants 
of collisions {i.e. t*(y), t*(4>(y; t*(y))), etc.) in the MSPD started at y, on the time interval 
[0,A]. For all r G {0, ...,R + 1}, let us denote y^ := 4>(y;tr) and y^ := 4>[A(y)](y; t^), so that 
yo = yo = y while yR+i = 4>A(y) and y^+i = 4>A(y)- Besides, the definition of the MSPD and 
flow property of both the TSPD and the MSPD yield, for all r G {1,..., R + 1}, 

Yr = 4>[A(y^_i)](y^_i;t^ - U-i), fr = 4>[A(y)](y^_i; - p-i). 


As a consequence. 


R -|-1 


II^A(y) - ^A(y)||i = Yi “ y’-lli “ “ y’-illi 

r=l 

R+l 

= E - ^[^(y)](yr-i;ir -ir-l)||l - ||yr-l -yr-l||l 

r=l 

1 R+1 


r=l 


i-keP4 


where the last line is obtained by applying the stability property of the SPD (1.2) typewise. Using 
the triangle inequality 


^fe(yr-l) - A^(y) < y(yr-l) - A^(yr-2) 


Afe(yi)-A^(yo) 


and summing by parts, we obtain 


II^A(y) -^A(y)||l < Y |A 2 (yr) - A;((yr-l) 

r=l 'y:keP4: 


Now the definition (1.8) of the coefficients combined with Assumption (LC) imply 


E |Afc(yr)-Ayyr-i) 

i-keP4 


r-f E 

j:keP4: '7^7 k' = l 


EJl/3:i(y)=^r}’ 


where {a : i,/3 : j) in the last term refers to (7 : : k') if 7 < 7 ' and to ( 7 ' : /c ',7 : k) if 

7 > 7 '. We note that the argument here is similar to the one employed in the proof of Fact 2 
in [13, Lemma 7.2.4]. 
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Since each pair [a \ \ j) G R(y) such that = U appears twice in the sum in the 

right-hand side above, we deduce that 


'y:k£P:^ k' = l (a:i,/3:j)GR(y) 




whence 


R 


Y1 Y1 |^fe(yd-A2(yr-i) 

r=l ^:keP^ 


< 


21/lc 


E 


(a:i,i3:j)eR(y) 


which completes the proof. 


□ 


For all y G we denote by 

NA(y) := y] l|^coii^^^(y)<A} 

(a:i,/3:i)eR(y) 

the number of collisions occurring in the MSPD started as y on the time interval [0, A]. Combining 
Lemma 4.2 with the first estimate (4.2), we obtain 

(4.4) ||$i(x)-$i(x)||i <2A:^^^NA($i-hx)), 

i=l 

for all L > 1. Of course, each term ^(^)) bounded by d{d — l)n^/2. Our purpose is to 

exhibit a bound of this order (with respect to n) for the whole sum, uniformly in L. We shall rely 
on the following remark, which is a consequence of Assumption (USH) and of the boundedness of 
the velocities. 

Remark 4.3. Let y G and a : p : j ^ such that a < p. Denote z := 4>A(y). 

• If (a : i, /3 : j) ^ R(y), then {a:i,P : j) ^ R(z). 

• d ~ - yj ~ Vi ~ .^ushA and, if (a ; ; j) e R(y) and T=°“^,j(y) < A, then 

o<ypyi < 2LcA. 

We deduce from the first part of the remark that, for all i G ,L}, R(4>^^(x)) C R(x), 

which allows us to rewrite 

L L 

^NA(4>i“^(x)) = E E l{(a:i,/3:j)eR(i.i-hx)),r-‘;^^.(i.i-hx))<A}- 

^=1 (ck::i,/3:jf)GR(x) £=1 

Now for all {a : i,/3 : j) G R(x), let £a:i,p:j be the lowest index in {1,... ,1/} such that {a : : 

j) G R(4>^^(x)) and < A. If there is no such index then the contribution of 

{a : i^/3 : j) in the sum above is null. Otherwise, the second part of the remark implies that after 
at most 

m := |"2 Lc/Lush1 

iterations, the pair (a : i,/3 : j) can no longer belong to R(4>^'"’^'^ ^^’^(x)). We deduce that 

L 

E^{(a:i,/3:i)eR(4i-^x)).T“"^..(4i-i(x))<A} ^ 

^=1 

and therefore conclude that 

L 

^Na(4>^^(x)) < md{d — l)n^/2. 
r=i 

Injecting this bound in (4.4), we complete the proof of (4.1) with the constant 
(4.6) C := d{d — 1)£iI/lc [^^c/^ushI • 
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Now, for t G [{L — 1 )A,I/A], by reasoning like in the derivation of (4.2), one obtains 

L-l 

||$(x;t) - $A(x;dl|i < E ll^t-^A($A - l>A)$i-'(x)||i 

£=1 

+ ll(^t-(L-l)A — ^t-(L-l)A)^A ^Wlll 

<AEll(^A-$A)i>i-'(x)||l 

1=1 

+ ll(^t-(L-l)A - 4’t-(L-l)A)4’A~^ Will- 

Since by Lemma 4.2, 

||($t_(L_l)A - ^t-(L-l)A)^i”^(x)||l < 2{t - (L - 1)A)^^ y] 

(a:i,/3:i)eR(y) 

(a:i,/3:i)eR(y) 

the above derived upper-bound for ||4>£(x) — 4>£(x)||i still holds for all t G [{L — 1)A,I/A]. 


4.2. Rate of convergence of the MSPD and its TSPD approximation with step A. 

The purpose of this subsection is to estimate the distance between the semigroup solution 
u = (r^, ..., u^) to the diagonal hyperbolic system (1.7) with initial data uq = (i7*m^,..., H^m^) 
for some m = (m^,... ,m^) G P(M)^, and the empirical CDF [x(n)] = «[x(n)],..., M^[x(n)]) 
associated with the MSPD started at some configuration x(n) G over finite time horizons. 

We also compare u with the empirical CDF Un,A[x(n)] = {u\ a [^(^)]5 • • • ?'^n,A[x(n)]) obtained 
with the iterated TSPD with step A starting from x(n), that is 

1 "" 

k=l 

Given x G we set /io[x] = (/iJM: • • • ^MoM) ^ P(I^)'^. 


Theorem 4.4. Under Assumptions (LC) and (USH), let n > 1, x(n) G and A > 0. Then, 
for all t > 0, 

||u(t, ■) - u„[x(n)](i, •)||Li(R)<i < Cl iJ,o[:>c{n)]) + ^ , 

||u(t, ■) - u„,A[x(n)](i, ■)||li(r)<j < £i ^wf^(m,/xo[x(n)]) + +CA, 

where the eonstant C is given in (4.6). 

When for a// 7 G {1,..., d}, m^([a^, b^]) = 1 with —oc < < +00 and 

(4.7) Vj-.keP^, xl{n) = {u2)~^ ( ^\n ’ 


then 


||u(i, •) - u„[x(n)](t, •)||li(r)<* < Cl 


ELi(&"'-«"')+^2diLC 


2n 


||u(t, •) - u„,A[x(n)](i, •)||Li(R)<i < C 


E7i(&"'- a"") + i2dLLC 


2n 


+ CA. 


Proof. The estimations involving u^^a are an immediate consequence of Proposition 4.1 and the 
estimations involving Un[x(n)]. To prove those estimations, following the approach described in 
the introduction of the article, we first write, for alH > 0 , 

\\u{t, ■) - u„[x(n)](i, •)||Li(R)<i < ||u(i, •) - Uoo[x(n)](i, •)||Li(R)<i 

+ ||uoo[x(n)](t,-) -u„[x(n)](i, •)||Li(R)<i, 


(4.8) 
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where Uoo [x(n)] refers to semigroup solution to the diagonal hyperbolic system with initial condition 
Uoo,o[x(n)] = {H * /io [^(^)]5 • • •, ^ * /^o[^(^)])- The stability property of Theorem 1.6 yields 

||u(i, ■) - Uoo[x(n)](i, ■)||li(k)‘* ^ -CiUuo - Uoo,o[x(n)]||Li(M)<i = £iWf^(m,/io[x(n)]), 

which, according to (hi) in Lemma 2.1, is smaller than ~ <^^)/(2n) if b'^]) = 1 for 

all 7 , and x(n) is given by (4.7). The following proposition allows to control the second term in 
the right-hand side of (4.8). □ 


Proposition 4.5. Under Assumptions (LC) and (USH), let n > 1 and x G D^. For all t >0, 

Proof. For all n > 1, let us extend the duplication operator introduced in the proof of Proposi¬ 
tion 3.2 by setting 

I Di ^ 

I X X 

where ^ 2 /c-i ~ ^ 2 k ~ all /c G {1,..., n} and 7 G {1,..., d}. Notice that 

(4.9) Vx,yGL>^, ||x-y||i = ||x-y||i. 


Let X G D^. Like in the proof of Proposition 3.2, we are going to estimate ||u 2 n[x](t, •) — 
Un[x](t, OIIl^r)^- The direct analysis of the MSPD being delicate, we introduce a ficticious step 
A > 0 to transform this analysis into the comparison between the TSPD with n particles and 2n 
duplicated particles on each time-step. Let t > 0 and L = [t/A]. One has 


L-l 


$(x;t) - i'(x;i) = ^(x) - W) + ^(x) - $t-(L-l)A^i ^(x). 


1=1 


Hence, by the triangle inequality and the stability of the MSPD, 

I|u2„[x](i,-) -U„[x](i,-)||L1(K)‘* = ll^(x;i) - ^(x;dlli 


L-l 


(4.10) 


<£i^||Mi-i(x)-$i(x) 




+ ll^t-(L-l)A^i ^(x) - i>t-(L-l)A^i "(x)||l. 


L-l/ 


Now, by the triangle inequality and (4.9) 

ii$A$Tkx) - ilwiii < ii^A^ykx) - $A$'fkx)iii 

+ ||l>A$f^)-II$i-'(x)||i 

+ |||.A$i-^(x)-$i(x)||i. 

The second term in the right-hand side is a comparison at time A between the TSPD with n 
particles starting from 4>^^(x) and the TSPD with 2n particles starting from the duplicated 

vector 4>^^(x). Reasoning like in the derivation of (3.3), we bound it from above by AdI/Lc/(2n). 
By Lemma 4.2, 

||$A$ykx) - $A$Tkx)||l + ||^A$i-'(x) - $i(x)||i 
< (^NA($Tkx)) +4NA($i-'(x))^ , 

where Na has been defined in (4.3). Plugging these estimations in (4.10) and dealing in the same 
way with the last term in the right-hand side of (4.10), we deduce that 

||u2„[x](t, •) - U„[x](i, •)||li(k)‘* ^ ^ (NA(^i”^(x)) + 4NA($i”^(x))) . 
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Clearly, for all a : j G with a < /3, 

L 

^=1 

We now use the same arguments as in the derivation of (4.5) to obtain that, for d\\ a \ P \ j G P^n 
with a < P, 


(4.11) 


L 




< |" 2 I/c/I/ush 1 - 


To this aim we only have to replace Remark 4.3 with the following: let y G and a : P : j ^ P^ 
such that a < /3^ and denote z := 4>A(y). 

• If (a : : j) ^ R(y), then for all (i',/) G {2i-l,2i} x {2j- l,2j}, {a : i',/3 : f) ^ R(z). 

• For all {i'J') G {2i - 1, 2i} x {2j - 1, 2j}, - z^, < - yf - Tush A and, if (a : i, ^ : 

j) G R(y) and < A, then 0 < y^ - yf < 2TcA. 

The remainder of the argument is the same and leads to (4.11). As a consequence. 


^(NA($i-'(x)) +4NA($i-'(x))) < 2d{d-l)n^{l + r2Lc/iusHl), 


£=1 


and 

||u 2 „[x](t, ■) - u„[x](t, ■)||Li(R)d < + Ad{d- 1)£iLlc( 1+ r2Lc/-LusHl)- 

Taking the limit A —> 0, we deduce that 

l|U 2 n[x](i,-) -U„[x](i,-)||L 1 (R)‘* ^ ' 

With this estimation replacing (3.3), we shall conclude like in the proof of Proposition 3.2. We 
therefore need to show that 


lim ||u2M„[x“](i, •) - Uoo[x](i, OllLifE)^! = 0. 

To this aim, we denote := /io[x] and recall the definition (1.9) of the discretisation operator to 
write 

||U2M„[x“](i, •) - Uoo[x](i, •)IIl 1(E)<* < l|U2M„[x“](i, •) - U2M„[x2M„m„](i, •)||l 1 (K)'‘ 

+ l|u2M„[X2Mnm„](t, •) - Uoo[x](t, OIIlHB)^!- 

Using the same arguments as in the proof of Proposition 3.2 but with Theorem 1.6 replacing 
Theorem 1.3, in particular the finite speed of propagation for both the MSPD and Uoo[x], we 
obtain that the second term in the right-hand side vanishes. On the other hand, the discrete 
stability estimate of Proposition 1.5 yields 

||u2M„[x^](i, •) - U2M„[x2M„m„](i, •)||Li(M)<i = ||$(x^;i) - $(x2M„m„; t)||i 

< Ullx^ - X2'^nm„||i 
= IIo[X2^n^n])- 

It is a property of the discretisation operator that, for all 7 G {1,... ,(i}, /iQ[x 2 ^ni^n] converges 
weakly to m'^ = /igM when M +00 [13, Lemma 8.1.5]. As a consequence, the corresponding 
empirical CDF converge dx-almost everywhere. Since the support of each measure /Tq [x 2 A^n^^n] 
is contained in the compact set [xj,x^], Lebesgue’s theorem implies the convergence in L^(M) of 
these CDFs, which implies 

lim w7(m„,//o[X2M„m„]) = 0 

M^ + OO 

by (1.10), and thereby completes the proof. □ 
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As an immediate corollary of Theorem 4.4, we obtain a convergence result for the MSPD to the 
semigroup solution of the system (1.7) which holds under less stringent conditions on the sequence 
of initial configurations. 

Corollary 4.6. Under Assumptions (LC) and (USH), let m G P(M)^ and let (x(n))n>i he such 
that for all n>l, x(n) G and 

(4.12) lim W^^\/ro[x(n)], m) = 0. 

n^+oo 

For all t > 0, the empirical CDF Un[x(n)](t, •) converges in L^(]R)^ to the semigroup solution u(t, •) 
of the system (1.7) with initial data uq = {H * ..., * m^). 

Under the condition (4.12), this corollary extends both [13, Theorem 2.4.5] and [13, Theo¬ 
rem 2.6.5]: in the former it allows to identify the limit, in the latter it relaxes the assumption that 
x(n) = Xn^- We however underline that a necessary condition for (4.12) to hold is that, for all 
7 ^ {1 ,..., d}, have a finite first-order moment. 

5. Numerical implementation 

5.1. Scalar conservation laws. This subsection is dedicated to the numerical implementation 
of the SPD in order to approximate the entropy solution to the scalar conservation law (1.1). The 
algorithm we use is described in §5.1.1. Then two case studies are presented: §5.1.2 addresses the 
Burgers equation 

(5.1) dtu + d,(^^=Q, 
with the CDF of the two-atom measure 

(5.2) m= i((5_i+(5i) 

as initial datum, while §5.1.3 deals with the conservation law with concave flux function 

(5.3) dtu + d, = 0, 

with the CDF of the two-sided exponential measure (or Laplace distribution) 

(5.4) m(dx) = - exp( —|x|)dx 
as initial datum. 

5.1.1. Numerical computation of the SPD. Given t > 0, a vector of (ranked) initial positions 
X = (xi, ..., x^) G Dn^ and a vector A = (Ai,..., A^) G of initial velocities, we use a remark 
due to Brenier and Grenier [6, Section 4] to devise an algorithm computing the vector 0[A](x;t) 
of positions at time t in 0{n) elementary operations, without following the detailed trajectory of 
each particle on the time interval [0 ,t]. 

Let us define the free transport flow '0[A](x; t) G by, for all /c G {1,..., n}, 

'0/c[A](x;t) := XkFtXk, 

and introduce the functions Ft and Qt on [0,1] such that Pt{0) = Qt(0) = 0 and, for all k G 
{l,...,n}, 

k k 

Pt{k/n) := (/)/c/[A](x;t), Qt{k/n) := Y 

k' = l k' = l 

with linear interpolation on [{k — l)/n,k/n]. Brenier and Grenier pointed out that Ft is the convex 
hull of Qt. Thus, our algorithm consists in computing the vector {Qt{k/n)^k = l,...,n} first, 
which obviously requires 0{n) operations, and then deducing Ft from the algorithm described by 
the pseudo-code below, which follows Andrew’s monotone chain algorithm [1], based on the Graham 
scan [10]. Note that the parallelisation of the computation of the vector {Qt{k/n), /c = 1 ,... ,n} 
is straightforward, while a parallelisable method to determine the convex hull of a sorted point set 
in the plane was devised in [9]. 
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The input is the vector Q of size n+1, with elements Q(k) = Qt{k/n) indexed by /c G {0,..., n}. 
The algorithm constructs a list L of integers k (ranked in decreasing order) containing the suc¬ 
cessive points at which Pt{k/n) = Qt{k/n). Once L is computed, Pt is reconstructed by linear 
interpolation, which finally yields 0[A](x;t). The length of the list L is denoted by |L| and its 
elements are indexed starting from 1. It is assumed that its first and second elements can be 
accessed and removed in constant time, 
initialise L = [1,0] 
for k = 2 to n 

while |L| > 1 and (Q(k)-Q(L(l)))/(k-L(l)) < (Q(L(1))-Q(L(2)))/(L(l)-L(2)) 
remove the first element from L 
end while 

insert k at the beginning of L 
end for 

It is easily checked, by induction on k G n}, that after the {k — l)-th iteration of the 

for loop, L contains the indices of the points at which Qt coincides with the convex hull of 
{Qt{k'/n), k' = 0,..., A:}. Indeed, the while loop consists in removing from L the points L(l) at 
which the piecewise linear function interpolating between the values of Qt at the points L(2), L(l) 
and k is concave, see Figure 1. 





Figure 1. An iteration of the algorithm computing L. The left-hand figure dis¬ 
plays the composition of L at the beginning of the (k-l)-th iteration. On the 
central figure, the point k-1 is removed from L. The right-hand figure displays the 
composition of L at the end of the k-th iteration. 

From the fact that the list L is browsed at each iteration, one could think that the algorithm uses 
0{in?) operations. This is actually not the case, as elements of L for which the test in the while 
loop returns true are removed from L and will therefore not appear again in the next iterations. 
As a consequence, the actual complexity of this algorithm is 0{n). 

Notice that an algorithm computing the explicit space-time points of collisions in the SPD 
would also take 0{n) operations, as there are at most n — 1 such points. As a consequence, the 
method presented here has the same computational efficiency as the explicit simulation of the SPD. 
However, the Brenier-Grenier trick allows for a significative simplification of the implementation. 

5.1.2. Burgers equation with two-atom initial measure. We consider the Burgers equation (5.1) 
with the CDF of the two-atom measure (5.2) as initial datum. In the SPD, two fans of particles 
are created, respectively originating from the points —1 and 1, see Figure 2 (a). These fans 
correspond to the fact that the entropy solution is the superposition of two rarefaction waves, 
respectively located at time t on [— 1 , — 1 -h t/ 2 ] and [1 - 1 - t/ 2,1 + t], see Figures 2 (b) and (c). 

The error between the particle approximation and the solution is plotted as a function of 
n, for several given terminal times t, on Figure 2 (d). In accordance with Proposition 3.2 and 
Remark 3.3, and since there is no discretisation error of the initial condition here (for even n), it 
is observed that this error is of the order of magnitude 1/n. 
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(a) SPD trajectory 
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Figure 2. Numerical results for the Burgers equation with two-atom initial mea¬ 
sure: (a) trajectories of 50 particles in the space-time plane, (b) value of the 
solution x) in the space-time plane, (c) profile of the solution u{t^ x) at succes¬ 
sive times t = 0, 8,16,..., 80 and (d) logarithmic plot of the error between 
the approximation obtained with 2^ particles and the solution, as a function 
of p = The different lines correspond to different values of t, namely 

t = 10, 20,..., 50, with the higher curves corresponding to the larger times. The 
slope of each line is —0.693 — log 2, which expresses the order 1/n of the error. 


5.1.3. Concave flux function. We now consider the conservation law with concave flux function (5.3) 
and the CDF of the two-sided exponential measure (5.4) as initial datum. As is made clear on 
Figure 3 (a), the particles progressively aggregate at 0. It results in the formation of a shock wave 
in the solution, see Figures 3 (b) and (c). The error is displayed on Figure 3 (d) and exhibits 
the following behaviour: given t > 0, there exists a critical number of particles such that: 

• below this number, the error does not vary with n, 

• above this number, the error decreases when n increases at the same rate as for the dis¬ 
cretisation of the initial measure. 

This behaviour is explained by the fact that, for n small, all the particles have arrived at 0 at time 
t, so that the approximate solution is the Heaviside function whatever n. But as soon as n is large 
enough to allow some particles to have an initial position far enough from 0 so that they have not 
reached 0 at time t yet, then the contribution of these particles in the approximate solution allows 
the latter to fit better the part of u outside of the shock wave, at the same rate as for the initial 
discretisation since the shape of u outside of the shock wave is merely an affine transformation of 
the initial profile. Following the conclusions of Section 2, this discretisation error is of the order 
log(n)/n. 
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Of course, the larger t is, the larger the uiagnitude of the shock wave is, therefore the better the 
Heaviside function approximates u and the more particles it takes to reach the critical number, 
which explains the ordering of the different curves on the picture. 



(c) Profile of u 



Position X 


(d) Error 



Figure 3. Numerical results for a concave ffux function with two-sided exponen¬ 
tial initial measure: (a) trajectories of 30 particles in the space-time plane, (b) 
value of the solution u{t^ x) in the space-time plane, (c) profile of the solution 
u(t^x) at successive times t = 0, 0.5,1,..., 5, (d) logarithmic plot of the er¬ 
ror between the approximation obtained with 2^ particles and the solution, as a 
function ofp= 1,...,12. The different lines correspond to different values of t, 
namely t = 0,l,...,17, and the higher curves correspond to the smaller times. 


5.2. Diagonal hyperbolic systems. We now turn to the numerical resolution of the diagonal 
hyperbolic system (1.7) thanks to the MSPD. 

A first method to simulate the trajectory of the MSPD obviously consists in computing the exact 
space-time position of each collision (between particles of the same type, or between particles 
of different types). The number of such collisions is at most of order n^: indeed, because of 
Assumption (USH), there are at most n^d{d — l)/2 collisions between particles of different types, 
and whenever particles of the same type collide, the space-time point of the next collision with a 
particle of another type is the same for all the particles in the current cluster. As a consequence, an 
algorithm computing the exact trajectory of each particle is expected to perform O(n^) elementary 
operations. We however believe that such an algorithm with optimal complexity would require a 
rather technical implementation. 

We therefore suggest to use a second method, which consists in approximating the MSPD by 
the iterated TSPD scheme described in Subsection 4.1. Thanks to the Brenier-Grenier algorithm 
introduced in the scalar case, each iteration of the SPD is easily implemented and requires 0{n) 










































Optimal convergence rate of the multitype sticky particle dynamics 


23 


elementary operations. Then we shall show below that updating the velocities after each step also 
requires 0{n) operations. As a consequence, computing the iterated TSPD on L steps requires 
0{nL) elementary operations. On the other hand, the error between the solution to the system (1.7) 
and the approximated solution provided by the iterated TSPD scheme was proved in Theorem 4.4 
to be of order 0(t/n +A) at time t. For the terms t/n and A to be of the same order, one therefore 
has to run the iterated TSPD scheme on L t/A n iterations, which leads to a total number 
of elementary operations in 0{n^). As a conclusion, this method has the same cost as the exact 
simulation of the MSPD, but it seems easier to implement. 

It follows from this discussion that to reach an approximation error of order e at time t > 0, the 
iterated TSPD scheme requires je^) elementary operations. In comparison, standard upwind 
schemes for the hyperbolic system (1.7), with a time step At and a mesh size Ax satisfying the CFL 
condition LqAI < Ax, are generally expected first-order accurate [14], so that the approximation 
error at time t writes C(t)Ax, with a constant C{t) that depends neither on At nor on Ax, and 
grows at least linearly with t. Besides, at each iteration of such a scheme, 0(1/Ax) elementary 
operations are necessary to compute the values of the fluxes and of the solution on the grid, so 
that after L tjAt iterations, 0{t/{AtAx)) elementary operations have been performed. As a 
consequence, the minimal number of elementary operations to reach a precision of order e at time 
t is obtained when the CFL condition is saturated, and it is in 0{tC{t)^ je^), which has the same 
dependence on e as the iterated TSPD scheme. 


5.2.1. Description of the algorithm. In order to simulate the MSPD, we use the approximation of 
the latter by the TSPD on small time steps A, as is described in Subsection 4.1. Given x G D^, we 
thus compute 4>^(x) instead of 4>(x;I/A). To this aim, we use an elementary iterative algorithm 
which will therefore perform L steps. At each iteration £ G {1,..., L}, it is necessary to: 

(i) compute the vector of velocities for the TSPD started from 4>^^(x), 

(ii) compute the evolution of each subsystem of particles according to the TSPD. 

Of course, the second step uses the algorithm described in Subsection 5.1 and therefore makes 
0{nd) elementary operations. The first step is realised by the following pseudo-code, the input 
of which is an array x of size d x n, such that x(gamma,k) contains the initial position xj. of the 
particle 7 : k. We recall that, for fixed 7 , x^ < for all /c G {1,..., n — 1}. 

current_indices = vector [0, 0] of size d 

while min(current_indices) < n 

gamma = max( argmin( x(g,current.indices(g)) ; 

g such that current.indices(g)<n ) ) 
k = current.indices(gamma) 
current.indices(gamma) = k+1 

velocity(gamma,k) = lambda(gamma,current.indices) 
end while 


In this pseudo-code, lambda (gamma, [k.l. 


jn 

n 

J W={kry — l)/n 



k_d]) returns the velocity 


^7—1 


n 


w. 


^7+1 

n 



dw 


of the particle j : k, so that at the end of the algorithm, the vector velocity (gamma, :) contains 
the initial velocities of the particles of type gamma for the SPD. 

There are nd iterations of the while loop, and to select gamma it is necessary to scan the vector 
cur rent, indices, which costs d operations. As a consequence, the computation of the velocities 
requires 0{nd?) operations. Since d is a physical parameter, we only consider the complexity with 
respect to the numerical parameter n and therefore conclude that the computation of 4>^(x) is 
made in 0{nL) operations. 


5.2.2. Case study: p-system. The p-system 


dfU = dxV, 

dtv X dx{p{u)) = 0 , 
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is a simple model for isentropic gas dynamics in one space dimension, where u is the specific volume 
of the gas and v is its velocity. The function p{u) determines the pressure in terms of the specific 
volume, and must generically satisfy p'{u) < 0 for all r > 0. In the sequel, we shall furthermore 
assume that there exists > 0 such that 


(5.5) 


[ ^-p'iu)du = 1 , 

Ju=0 


which is the appropriate condition to develop our probabilistic approach, and in which case the 
specific volume will take its values in [ 0 , z^]. 

Let us define the Riemann invariants w~ and w~^ for this 2 x 2 system [15] by 

= V ± g(u)^ 

where, for all u G [ 0 , z/]. 


ru 2 

g{u):= c(r)dr--, c{u) := ^J-p’{u). 

Jr=0 ^ 


The assumption (5.5) ensures that, if w G [0,1], then {w~^ — w )/2 belongs to the image of 
g, and it is immediately checked that u and v are recovered from the formulas 

r)+ . 


(5.6) 


u = g 


-1 




w 


V = 


Furthermore, the Riemann invariants satisfy dfW^ = :^c{u)dxW^, which rewrites under the form 
of the 2 x 2 diagonal system 

( dtw~ + \~{w~^w^)dxW~ = 0, 

1 dtw^ + {w~ ,w^)dxw^ = 0 , 


(5.7) 

with 


A^(r; ,'zn+) :=^c(g 


-1 


= T 


(g-iy((nj+-w-)/2)' 


Under the assumption that p' be continuous and negative on [0, z/], one can define £ G (0, Too) by 

£ := inf c(r), 

0<u<u 

and get 

(5.8) Vre”, G [0,1], {w~^w~^) <—£<£< X~{w~ 

so that, for initial conditions Wq and Wq given by the CDFs of probability measures, the system 
satisfies Assumption (USH) with constant Tush = 2^. 

We now present numerical approximations of u and v for the choice of pressure function 


p{u) = 


4z^argsinm(zi:/2) 


arctan ( — 1 — arctan ( ZD ( —- 


(f) 


u e [o,v], 


where z^ > 0 is a given reference specific volume and /^ > 0 is a dimensionless shape parameter. 
This choice implies 

argsinh (k - |)) 


9{u) = 


2 argsinh 


SO that 


X^{w , w^) = ^ 


2 z/argsinh(zi:/ 2 ) cosh [{w~^ — w~) argsinh (|)] 


The relation (5.6) yields 


u = V \ - H— sinh {w^ — w ) argsinh j , 


We first address the case where the initial conditions Wq and Wq are the respective CDFs of 
the shifted two-sided exponential distributions m~ and defined by 


rri^{dx) := - exp(-|a; ± a;o|), 
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for some xq > 0. In this case, Wq{x) > Wq{x) for all x G M, which implies that •) — w~{t, •) 

remains nonnegative at all times t > 0. This is indeed easily checked at the level of the MSPD, a 
trajectory of which is plotted on Figure 4: if Wq and Wq are the empirical CDFs of two vectors 
= {xf^... G Dn^ then Wq > Wq if and only if, for all /c G {1,... ,n}, x^ < x^. By (5.8), 
for alH > 0 we have 

$+ (x; t)<x+ - £t, (x; t) > + it, 

with the obvious notation x = (x“,x+) G so that the corresponding empirical CDF satisfies 
•) > w~{t, •). That this inequality still holds true in the limit n +oc can be checked using 
the notion of trajectories introduced in [13, Section 5], see in particular [13, Corollary 5.1.2]. 

One can observe on Figure 4 that particles of the same type never collide with each other. This 
is due to the fact that, for fixed w~ G [0,1], the mapping i-^ A+(re“,re+) is increasing on 
[re“, 1 ]. As a consequence, two consecutive particles of type + with no particle of type — between 
them can only have velocities taking them away from each other. The same phenomenon occurs for 
particles of type —. However, collisions between particles of different types modify the velocities 
of these particles. Thus, the Riemann invariants w~ and respectively plotted on Figures 5 
(a) and (b), undergo two interacting rarefaction waves, drifting away from each other on account 
of Assumption (USH), without forming any shock. The specific volume u and the velocity v are 
finally plotted on Figures 5 (c) and (d). 



Position X 


Figure 4. Trajectory of the MSPD (obtained with the iterated TSPD scheme 
with A = 0.03) with 20 particles per type associated with the p-system for Wq > 

Wq . Blue rays correspond to particles of type —, red rays correspond to particles 
of type +. Here xq =0.1 and u = 0.5, k = b. 

As a sticky particle dynamics where particles never stick together may seem a little disappoint¬ 
ing, we now choose initial conditions that do not satisfy the condition that Wq (x) > Wq (x) for all 
X G M. To this aim, we still assume Wq to be given by the CDF of m~{dx) = exp(—|x — xo\)/2 
with Xq < 0, and take [x) = l| 2 ,>o}. Particles of both type can now aggregate into clusters, as 
is depicted on Figure 6 . But on account of Assumption (USH), after a finite time (that generally 
depends on the number of particles), the property that > w~ is recovered and the particles 
start drifting away from each other again. This is also observed on Figure 6 , where two clusters 
blow up under the effect of a collision. As a result, the functions w~ ^ u and v exhibit shocks on 
short times, and are essentially given by interacting rarefaction waves on longer times, see Figure 7. 
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(a) Riemann invariant w- (b) Riemann invariant w+ 



Position X Position x 


Figure 5. Representation in the space-time plane of the Riemann invariants w~ 
(a) and (b), of the specific volume u (c) and of the velocity v (d) for the 
p-system with Wq > Wq . The simulation uses the iterated TSPD algorithm with 
parameters n = 200 and A = 0.03, and u = 0.5, 6: = 5. 


Appendix A. Proof of Proposition 1.4 


Since (t,x) i-^ u{t + s^x) (resp. (t,x) i-^ v{t + s,x)) is the entropy solution to (1.1) with initial 
condition x u{s,x) (resp. x v(s,x)), it is enough to deal with the case 5 = 0. 

We define m and m' in P(M) by rq = ^ '^o = ^ aiid use the discretisation of m and 

m' corresponding to (1.9), namely 

/*(2i+l)/(2(n+l)) />(2i+l)/(2(n+l)) 

Xi(n) = (n-h 1) / UQ^(w)dw, p^(n) = (n + l) / VQ^(w)dw, 

Jw=(2i-l)/(2(n+l)) Jw=(2i-l)/(2(n+l)) 

for all i G {l,...,n}. Let x(n) = (xi(n),..., Xn(n)) and y(n) = (pi(n),..., p^(n)). According 
to [13, Lemma 8.1.5], /io[x(n)] (resp. jiio[y{n)]) converges weakly to m (resp. m') as n ^ +oo. 
Moreover, by [13, Lemma 8.1.6], 


By (1.2), for 


lim Wi(/io[x(n)],/io[y(^)]) = Wi(m,m'). 

n^+cxD 



and 
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Figure 6. Trajectory of the MSPD (obtained with the iterated TSPD scheme 
with A = 0.02) with 20 particles per type associated with the p-system for initial 
conditions with shocks. Blue rays correspond to particles of type —, red rays 
correspond to particles of type +. Red particles first remain aggregated into a 
single cluster up to the collision with the median blue particle, which makes it 
blow up then. Here xq = —4 and u = 0.5, n = b. 


one has 


\\(p[X]{x{n);t) - < Wi{iJ.o[x{n)], iJ.o[y{n)]) + t 

< Wi{po[x{n)], iio[y{n)]) + t 


'f' pi/n 

“ Kw))<iw 

/ \X{w) — jj.{w)\dw. 

J w=0 


One concludes by taking the limit n +oc in this inequality since Theorem 1.3 and the lower 
semi-continuity of Wi with respect to the weak convergence topology [17, Remark 6.12] ensure 
that 


||u(i, •) -f(t, OIIl^R) ^ liminf ||^[A](x(n);i) - (j)[iJ,]{y{n);t)\\i. 
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Figure 7. Representation in the space-time plane of the Riemann invariants w~ 
(a) and (b), of the specific volume u (c) and of the velocity v (d) for initial 
conditions with shocks. The simulation uses the iterated TSPD algorithm with 
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